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A new approach to the formulation 
of scalar flux closure 

By Y. Shabany AND P. Durbin 
1. Motivation and objectives 

The solution of fluid dynamics equations for a turbulent flow requires the mod- 
eling of turbulence statistics if the averaged form of the equations are used. This is 
usually the case except in direct numerical solution methods which are limited to 
low Reynolds numbers. The major effort of the researchers in this field is to develop 
closure models that improve the accuracy of turbulent flow predictions. However, 
it is understood that the more accurate the models are, the more complex they 
will be. The second order closure models seemed to provide a compromise between 
complexity and accuracy. In this class of models the exact equations for Reynolds 
stresses and scalar fluxes are derived and the unknown terms are modeled in terms 
of the other known parameters. 

The modeling of Reynolds stress and scalar flux transport equations is done sepa- 
rately, although the same approaches are used in most cases. It must be mentioned 
that the area of scalar transport has received less attention than momentum trans- 
port (Reynolds stress). Therefore, turbulence models for scalar fluxes are rather 
less well developed than models for Reynolds stresses. This may be in part because 
prediction of the mean flow and Reynolds stress is often a prerequisite to prediction 
of convective scalar transport. But, conversely, because of this intimate coupling 
between momentum and scalar flux, models of scalar transport may provide con- 
straints on the momentum model. 

This report shows that if a stochastic differential equation (Langevin equation) 
for velocity fluctuation vector is known, it is possible to derive the equations for 
scalar flux transport. Durbin and Speziale (1994) showed that the second moment of 
this stochastic differential equation gives an equation for the evolution of Reynolds 
stress tensor. Similarly, the stochastic equation will give an equation for scalar flux. 
Therefore, a coupling between these two is present. The basis for the present work 
is that there should be Langevin equations that can produce acceptable models for 
both the Reynolds stress tensor and the scalar flux vector. Having found this basic 
Langevin equation, the amount of work needed to model the second order closure 
problems is reduced; using the well developed models for Reynolds stress equations, 
it will be possible to derive corresponding models for scalar flux equation. 

2. Accomplishments 

2.1 Langevin equation and scalar flux closure 
The simplest Langevin equation for a random velocity vector is 

du % = -Uidt + y/c^edWi(t), 


( 1 ) 
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where u, is the velocity fluctuation vector, t is the independent variable time, T 
is the turbulent time scale ( k/e ) where k is the turbulent kinetic energy per unit 
mass and t is the rate of dissipation of k, and W)(f) is the Wiener stochastic process 
(Arnolds, 1974). ci and c 0 are constants which are determined later. It was shown 
by Durbin and Speziale (1994) that the second moment of this equation, 


duTuj 

dt 


+ coetijj, 


( 2 ) 


is an equation for the evolution of Reynolds stresses in the absence of mean velocity 
gradient provided that c\ = 2 and Co = 2/3. 

In homogeneous turbulent flow the position of a fluid particle is determined by 
the following equation. 

^.=u i +X j U iJ (3) 

Here, X{ is the Lagrangian position vector of the particle and Ui is the velocity of 
the fluid particle, which is at position X x at time t; U x j is constant in homogeneous 
turbulence. The dispersion tensor is defined as 


K xj =u i X J . 


(4) 


It can be shown that if the molecular diffusion of the scalar contaminant 0 is 
neglected (high Peclet number), the turbulent scalar flux is related to the dispersion 
tensor by _ ___ 

Ui0 = — (5) 

Therefore, if a transport equation for K XJ is known, the equation for the transport 
of scalar contaminant can be derived using Eq. 5. 

The transport equation for K x j is simply obtained by substituting Eqs. 1 and 3 
in 

d(u t Xj) = (u, + du x )(Xj + dXj) - u x Xj ^ 

= u x dXj + Xjdui + du x dXj, 

and averaging. This is the same method used by Durbin and Speziale (1994) to 
derive the transport equation for Reynolds stresses. The result is: 


dK t j ci 

- 2 T K “- 


(7) 


Note that the mean velocity gradient and therefore the second term of Eq. 3 is zero 
for the case considered here. The coefficient of the second term of Eq. 7 (ci/2) does 
not agree with the empirical values which are about 2cj where cj = 1.8 (Launder, 
1978). Therefore, Eq. 1 can not be used as a base Langevin equation for both 
Reynolds stress and scalar flux closure models. However, a modified form of this 
equation given as 


du, = 


-^-Uidt + \J ( 2c*v? — 1 )tj>i k dW', k + ^/cofdWt, 


( 8 ) 
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provides a consistent Langevin equation for both momentum and scalar flux trans- 
port. In Eq. 8 pij is the generalized square root of b %] defined by p] - ~p 2 kk 6 tJ = b % 

and c 0 = 3 [2cm — 1 — (cm — )p\k\ • Note that W t and W t ' are independent 

Wiener processes and dW t dW t = 0. It can be shown that the second moment of 
Eq. 8 is Eq. 2 and the evolution equation for Kij is 

d&ij cm 

-dT = -jrK ij + u i u j (9) 

The scalar gradient evolves by d{Q tj )/dt = -U k>J Q k (which is zerojn this case 
but not in general). Hence the transport equation of the scalar flux, u t 8, is simply 
obtained by taking d/dt of Eq. 5, substituting Eq. 9 for dK XJ /dt, and using the 
above mentioned evolution equation for 0 i; . The final result is 


* _ ^ M 7T 

= -UiUjQj - —u,9. (10) 

This is the equation for the transport of scalar flux in the absence of mean velocity 
gradient. The first and second terms are the production by mean scalar gradient 
and the slow part of the pressure-scalar gradient correlation respectively. 

The importance of this method is that there is no need to develop a separate 
closure model for the equation of scalar transport if there is already a closure model 
for the transport of Reynolds stresses. It was shown (Durbin and Speziale, 1994) 
that for any Reynolds stress closure model there is a Langevin equation, the second 
moment of which is that model equation. Having this Langevin equation, it is 
possible to derive a transport equation for the scalar flux by the method outlined 
above. 


2.2 Results 

The general linear model for the evolution of Reynolds stress tensor is 


du t Uj 

dt 


~ j( u ‘ u J - - c 2 (Pij - | PSij ) - c 3 (D tJ - 2 -P6 tj ) 

— c s kS t] + P ,j — -e6i j 


( 11 ) 


where 

Pij = —Ui U k Uj yk — UjU k Ui ,Jt? 
Dij = -UiUkU k j - uJu^U kti , 



A special case of this model is the IP model where c 3 = C s = 0. It can be shown 
that the second moment of Langevin equation 


du ‘ = -jruidt + y/{2 cm - crfepitdW;* + (c 2 - 1 )u k d k U,dt + y/E^dW t , (12) 
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is the special case of Eq. 11 corresponding to the IP model, provided that 


Co = 3 


2cm - 1 - c 2 y - (cm - 


Following the procedure mentioned in Section 2, it can be shown that the evolution 
equation for K{j is 


^ - -M-Kij + KikUj'k + (C2 - l)K kj U itk . 

at 1 


(13) 


Note that K tJ is not a symmetric tensor. The scalar flux equation is 


= ~^-u t 6 - UiUj&j + (c 2 - l)u k 9Ui y k- 

at 1 


(14) 


The significance of this result is that the coefficients of the scalar flux model are 
not independent of those in the Reynolds stress equations (Eq. 11). 

The dimensionless dispersion tensor is defined as 

D *j ~ k2 Kij ' 

An evolution equation for Dij can be obtained. However, in equilibrium, the rate 
of change of D t j is zero. Therefore, the following algebraic equation is obtained. 


D 


*; _ 


9k 


= Tij + DikiSlj-ulj) + (c 2 - 1 )(S* k +u,* k )D kj 


(15) 


where rij = u t uj/k , S*j — TSi ; , = Tu>,j and 


9k 


P P 

CM + c ( 2 — C (1 — |- 2 2 

€ € 


-1 


In a two-dimensional uniform shear flow: 

s ''=K° o) md “"=K-°s °) 

It can be shown that in this case 

D\\ = g k [^n + C 2 g k S*T \2 — 2(1 — c 2 )(0fcS*) 2 T2 2 ] , 
D \2 = g k [ri 2 + (c 2 - 1)0*S*T 22 ] , 

D 2 I = 9k [ r 12 + 9 k T 22 ] , 

D 22 = 9kT 22 * 


(14) 
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FIGURE 1. Eddy diffusivity profiles according to the IP model; , model 

prediction; *, DNS(HS) Pr=0.71; +, DNS(HS) Pr=2.0; o, DNS(HW) Pr=0.71; x, 
DNS(HW) Pr— 2.0. 


Evaluation of D 22 with cm — 3.4 is compared to the numerical data of Kim and 
Moin (1989) in Fig. 1. The DNS data are for heat source (HS) and heated wall (HW) 
cases. With cm = 3.4 the value calculated for D \ 2 is not in good agreement with 
the numerical data; a value close to 0.85 gives more reasonable results. However, 
this value of cm is preferred in order to predict the transverse scalar flux, v0, as 
accurate as possible. This component of scalar flux has the main contribution in 
channel flow and boundary layer heat transfer. 

The same calculations were done for the general linear model and the results are 


D 1 1 
D12 


9_k_ 

A 

9k 


r n + + c 2 — 2 )gkS*ri2 - 2(1 — c 2) ^*^ ^ r 22 


-(1 -C 2 )gkS*T22], 


D21 — 


9k 


D 22 = ^ [r 2 2 + c 3 g k S*Ti 2 ] , 


r l2 + C3gkS*Tn + (“ — l)gkS*T22 + 2C3 — i ^ T\ 2 


(15) 


where A = 1 + (1 - c 2 )c 3 (g k S*) 2 and g k = [c M + (2 - c ( \) f + c f2 - 2] \ The 
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expression for D 22 given by Eq. 15 shows that for high enough values of S* there is 
a possibility for D 22 to become negative if C3 is not equal to zero; for the channel 
flow case this happens in the near wall region. The negative value of D 22 does 
not have any physical interpretation and is not supported by DNS data of Kim and 
Moin (1989). However, a non- zero value of C3 is necessary to predict different values 
for 622 and 633 in homogeneous shear flow. 

4. Future work 

The main purpose of this research is to obtain a consistent way of deriving both 
Reynolds stress tensor and scalar flux vector closures from same Langevin equation 
for velocity fluctuation vector. The following main problems must be resolved before 
this goal is achieved: 

1. The coefficient of the slow term (cm) in the evolution equation of different com- 
ponents of scalar flux must be different in order to get a good agreement with 
experimental or DNS data for all the components of scalar flux. Therefore, a 
simple constant value does not seem to solve the problem. 

2. As mentioned at the end of Section 3 a non- zero value of C3 is necessary to 
differentiate the values of 622 and 633 in a homogeneous shear flow. On the other 
hand, a non-zero value of this constant causes the model to predict negative values 
of D 22 for high enough S* which seems unreasonable. 

The solution of these two problems is the main focus of this research. 
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